
###### INITIALIZE
setwd("C:/Users/Philipp/Projects/geomatch")
source("R/PGHandlerFunctions_20131001.R")
source("R/GeoDBMover_v1.0.R")
source("R/PGHandlerFunctions_20131001.R")
library(mgcv)
library(Hmisc)
library(reshape2)
library(spdep)
library(MASS)


###### FUNCTIONS
deleteOverlaps <- function(input.df) {
  # Get points without overlaps
  single.points <- unique(input.df$rpid[is.na(input.df$rpid2)])
  
  # Iterate through points with overlaps and randomly eliminate overlapping points (i.e., buffers)
  input.op.df <- input.df[!is.na(input.df$rpid2), names(input.df) %in% c("rpid", "rpid2")]
  # Shuffle data frame
  input.op.df <- input.op.df[sample(1:nrow(input.op.df)),]
  # Iterate through overlapping points and eliminate overlaps
  op.cand.points <- unique(input.op.df$rpid)
  op.keep.points <- op.cand.points
  for (i in 1:length(op.cand.points)){
    thispoint <- op.cand.points[i]
    if (thispoint %in% op.keep.points){
      oppoints <- input.op.df$rpid2[input.op.df$rpid == thispoint]
      print(length(oppoints))
      flush.console()
      op.keep.points <- op.keep.points[!(op.keep.points %in% oppoints)]
    }
  }
  
  return(c(op.keep.points, single.points))
}



###### CONNECT TO DB
con <- getPGConn("growup", 5432, "cederman.ethz.ch", "admin", "hNo7Yoo")



###### GET AND MODIFY BUFFER SAMPLE

# Get data from DB
rb.df <- dbGetTable(con, "hunzikp", "randombuffersdata")

# Data manipulation
rb.df$lnpop <- log(rb.df$pop1990 + 1)
rb.df$lnarea <- log(rb.df$area_sqkm)
rb.df$lncap <- log(rb.df$capdist_km + 1)
rb.df$lnborder <- log(rb.df$border_km + 1)
rb.df$cowid <- as.factor(rb.df$countries_cowid)
rb.df$groupcount09[is.na(rb.df$groupcount09)] <- 0
rb.df$groupcount65[is.na(rb.df$groupcount65)] <- 0
rb.df$gpdiff <- rb.df$groupcount09 - rb.df$groupcount65
rb.df$mgs09 <- ifelse(is.na(rb.df$min_groupsize09), 1, rb.df$min_groupsize09)
rb.df$max_noauton09 <- ifelse(rb.df$max_discrim09 == 1 | rb.df$max_powerless09==1, 1, 0)
rb.df$maxnoauton09[is.na(rb.df$max_noauton09)] <- 0
rb.df$max_excl09[is.na(rb.df$max_excl09)] <- 0
rb.df$max_excl[is.na(rb.df$max_excl)] <- 0
rb.df$max_powerless[is.na(rb.df$max_powerless)] <- 0
rb.df$max_discrim[is.na(rb.df$max_discrim)] <- 0
rb.df$mno <- ifelse(rb.df$max_powerless + rb.df$max_discrim > 0, 1, 0)
rb.df$lec <- log(rb.df$ethnologue_count+1)

# Subsetting
rb.df <- subset(rb.df, (ssafrica == 1))  # Regional subsetting
rb.df <- subset(rb.df, cntr_petropoint == 1)  # Only countries with petroleum fields
rb.df <- rb.df[!is.na(rb.df$lnpop),]
rb.df <- rb.df[!is.na(rb.df$lnborder),]


###### RANDOM THINNING WITHOUT MATCHING

# Get table with overlaps for petropoints
op.df <- dbRunScript(con, "SQL/getBufferOverlaps_v0.1.sql", return=TRUE)

# Get IDs for randomly thinned sample
no.points <- deleteOverlaps(op.df)

# Get randomly thinned data frame
th.df <- rb.df[rb.df$rpid %in% no.points,]



###### UPLOAD MATCHED DATA AND GET WEIGHTS

retval <- dbCreateTable(con, "hunzikp", "rbmatch", "rpid", "int", drop=TRUE) 
retval <- dbInsertDF(con, "hunzikp", "rbmatch", th.df[,"rpid",drop=FALSE])

weights.df <- dbRunScript(con, "SQL/getWeights_v0.1.sql", return=TRUE)
th.df <- th.df[order(th.df$rpid),]


# Within country all weights matrix
cw.df <- subset(weights.df, select=c("rpid", "rpid2", "cdist_km"))

cw.wide.df <- dcast(cw.df, rpid ~ rpid2)
cw.mat <- as.matrix(cw.wide.df[,-1])
rownames(cw.mat) <- colnames(cw.mat) <- unique(cw.df$rpid)

cw.list <- mat2listw(cw.mat, row.names = rownames(cw.mat), style="W")
cw.nb <- cw.list$neighbours
cw.wmat <- nb2mat(cw.nb, zero.policy=TRUE)
cw.smat <- Matrix(cw.wmat, sparse=TRUE)


# Within country top 3 weights matrix
t3cw.df <- subset(weights.df, select=c("rpid", "rpid2", "top3cdist_km"))

t3cw.wide.df <- dcast(t3cw.df, rpid ~ rpid2)
t3cw.mat <- as.matrix(t3cw.wide.df[,-1])
rownames(t3cw.mat) <- colnames(t3cw.mat) <- unique(t3cw.df$rpid)

t3cw.list <- mat2listw(t3cw.mat, row.names = rownames(t3cw.mat), style="W")
t3cw.nb <- t3cw.list$neighbours
t3cw.wmat <- nb2mat(t3cw.nb, zero.policy=TRUE)
t3cw.smat <- Matrix(t3cw.wmat, sparse=TRUE)


# Across countries top 3 weights matrix
t3w.df <- subset(weights.df, select=c("rpid", "rpid2", "top3dist_km"))

t3w.wide.df <- dcast(t3w.df, rpid ~ rpid2)
t3w.mat <- as.matrix(t3w.wide.df[,-1])
rownames(t3cw.mat) <- colnames(t3w.mat) <- unique(t3w.df$rpid)

t3w.list <- mat2listw(t3w.mat, row.names = rownames(t3w.mat), style="W")
t3w.nb <- t3w.list$neighbours
t3w.wmat <- nb2mat(t3w.nb, zero.policy=TRUE)
t3w.smat <- Matrix(t3w.wmat, sparse=TRUE)


# Across countries all weights matrix
w.df <- subset(weights.df, select=c("rpid", "rpid2", "dist_km"))

w.wide.df <- dcast(w.df, rpid ~ rpid2)
w.mat <- as.matrix(w.wide.df[,-1])
rownames(w.mat) <- colnames(w.mat) <- unique(w.df$rpid)

w.list <- mat2listw(w.mat, row.names = rownames(w.mat), style="W")
w.nb <- w.list$neighbours
w.wmat <- nb2mat(w.nb, zero.policy=TRUE)
w.smat <- Matrix(w.wmat, sparse=TRUE)




#### GROUP COUNT

## IID COUNT MODELS
# DV: groupcount09
qpois.rt.fit <- glm(groupcount09 ~ petropoint, data=th.df, family="quasipoisson")
summary(qpois.rt.fit)
qpois.cov.rt.fit <- glm(groupcount09 ~ petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec + cowid, data=th.df, family="quasipoisson")
summary(qpois.cov.rt.fit)


## SPATIAL REGRESSION
# DV: groupcount09
slm.t3cw.fit <- lagsarlm(groupcount09 ~ petropoint, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- lagsarlm(groupcount09 ~ petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- lagsarlm(groupcount09 ~ cowid + petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- errorsarlm(groupcount09 ~ cowid + petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)



#### EXCLUSION

## IID MODELS
# DV: Mean Exclusion
lm.fit <- glm(mean_excl ~ petropoint*lnborder + groupcount09*lnborder, data=th.df)
summary(lm.fit)
lm.cov.fit <- glm(mean_excl ~ petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec + groupcount09, data=th.df)
summary(lm.cov.fit)
lm.cov.fit <- glm(mean_excl ~ cowid + petropoint*lnborder + lnarea + lnpop + lncap + lnborder + elevsd + lec + groupcount09, data=th.df)
summary(lm.cov.fit)

# DV: Mean No Auton
lm.fit <- glm(mean_noauton ~ petropoint + groupcount09, data=th.df)
summary(lm.fit)
lm.cov.fit <- glm(mean_noauton ~ petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec + groupcount09, data=th.df)
summary(lm.cov.fit)
lm.cov.fit <- glm(mean_noauton ~ cowid + petropoint*lnborder + lnarea + lnpop + lncap + lnborder + elevsd + lec + groupcount09, data=th.df)
summary(lm.cov.fit)



## SPATIAL REGRESSION
# DV: Exclusion 09
slm.t3w.fit <- semprobit(max_excl09 ~ petropoint, data=th.df, W=t3w.smat, ndraw=100, burn.in=10, showProgress=TRUE)
summary(slm.t3w.fit)
slm.t3w.fit <- semprobit(max_excl09 ~ cowid + petropoint + log(ethnologue_count+1) + lnarea + lnpop + lncap + lnborder + elevsd, data=th.df, W=Matrix(t3w.wmat, sparse=TRUE), ndraw=500, burn.in=100, showProgress=TRUE)
summary(slm.t3w.fit)

# DV: Mean Exclusion
slm.t3cw.fit <- errorsarlm(mean_excl ~ petropoint*lnborder + groupcount09*lnborder, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- errorsarlm(mean_excl ~ petropoint*lnborder + lnarea + lnpop + lncap + lnborder + elevsd + lec + groupcount09*lnborder, data=th.df, listw=cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- errorsarlm(mean_excl ~ cowid + petropoint*lnborder + lnarea + lnpop + lncap + lnborder + elevsd + lec + lnborder*groupcount09, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)

# DV: Mean No Auton
slm.t3cw.fit <- errorsarlm(mean_noauton ~ petropoint + groupcount09, data=th.df, listw=cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- errorsarlm(mean_noauton ~ petropoint + lnarea + lnpop + lncap + lnborder + elevsd + lec, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)
slm.t3cw.fit <- errorsarlm(mean_noauton ~ cowid + petropoint*lnborder + lnarea + lnpop + lncap + lnborder + elevsd + lec + groupcount09*lnborder, data=th.df, listw=t3cw.list, zero.policy=TRUE)
summary(slm.t3cw.fit)



